Fractal dimension of domain walls in the Edwards- Anderson spin glass model 
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We study directly the length of the domain walls (DW) obtained by comparing the ground states of the 
Edwards- Anderson spin glass model subject to periodic and antiperiodic boundary conditions. For the bimodal 
and Gaussian bond distributions, we have isolated the DW and have calculated directly its fractal dimension d t. 
Our results show that, even though in three dimensions df is the same for both distributions of bonds, this is 
clearly not the case for two-dimensional (2D) systems. In addition, contrary to what happens in the case of the 
2D Edwards-Anderson spin glass with Gaussian distribution of bonds, we find no evidence that the DW for the 
bimodal distribution of bonds can be described as a Schramm-Loewner evolution processes. 
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PACS numbers: 75.10.Nr, 75.40.Mg 

Whereas the properties of long range Ising spin glasses are 
now well understood, after more than 20 years of work, the 
same cannot be said about short range spin glasses. This 
is true even for very simple models: the nature of the or- 
dering of the low temperature phase of the two-dimensional 
(2D) Edwards- Anderson (EA) spin glass model 1 is still be- 
ing debated. Even though the fact that at T = EA mod- 
els with Gaussian (EAG) and bimodal (EAB) bond distribu- 
tions belong to two different universality classes seems well 
established, 2 new studies 3 - of low energy excitations (fractal 
droplets) show that there is still room for discussion. 

It has recently been suggested 4 ' 5 that domain walls (DWs) 
can be described as Schramm-Loewner evolution (SLE) pro- 
cesses. These processes are Brownian walks of diffusion con- 
stant k and fractal dimension df = 1 + fc/8. Furthermore, 
using conformal field theory the stiffness exponent 9 (which 
characterizes the scaling of the DW energy) can be related to 
the fractal dimension d f via 



df E = 1 



(1) 



4(3 + 6)' 

This seems to be true for the 2D EAG, as for 9 = -0.287 (4 )& 
Eq. 0} gives d f SLE = 1.2764(4), which is compatible with 
the best numerical estimate df = 1.274(2). 7 It is not clear, 
however, whether such a relation should hold for the 2D EAB, 
because of the high degeneracy of its ground state (GS). If it 
did, using the fact that the stiffness exponent seems to vanish, 8 
Eq. Q) would yield d f SLE = 1 .25. 

In the EAB model the degeneracy of the GS precludes a 
clear-cut definition of the fractal dimension of the DW. For 
this reason, most of the estimates of df are based on a scaling 
argument of Fisher and Huse, 9 which states that the entropy 
of droplets of size L should scale as Sow ~ L d fl 2 . It must 
be stressed that this was originally proposed for systems with 
only one GS. The estimates obtained using this scaling range 
from d f w 1 .0 to d f = 1 .30(3)^11 Even though very recently 
more direct measurements have been attempted , 12 ' 13 in those 
works the sampling of the DWs was not controlled. 7 In Ref. 7 
this problem is avoided and bounds are provided for the true 
df. 1.095(2) <d f < 1.395(3). 

In this paper we present the results of an extensive numeri- 



cal study of the fractal dimension of DWs, using a direct mea- 
sure of their length. We have studied both 2D and 3D sys- 
tems with Gaussian and bimodal distributions of bonds, but 
we have concentrated on the EAB model. For small systems 
we have calculated the exact average DW length. For larger 
systems an estimate of this quantity has been calculated by 
choosing only one pair of GSs. Even though the sampling we 
use is clearly not uniform, the values we obtain with this es- 
timate coincide with the exact ones, for small systems. Our 
results show that, whereas in 3D df is the same for EAB and 
EAG, in the 2D case the corresponding values are clearly dif- 
ferent and, in particular, for the EAB, different from the value 
that would be obtained from Eq. (Q]i (assuming that 9 — 0). In 
addition, we have performed one test to see whether DWs in 
2D can be described as SLE processes. Even though for the 
EAG the result of this test is positive, for the EAB the outcome 
of the test is clearly negative. 

We start by considering the Hamiltonian of the EA 
model for spin glasses 1 on square and cubic lattices, H = 
Tt(ij)Jij a i°j> wnere Oi = ±1 is the spin variable and 
indicates a sum over nearest neighbors. The coupling con- 
stants are independent random variables chosen either from a 
bimodal (±7) or Gaussian bond distributions, both with zero 
mean and variance one. In 3D we have used periodic bound- 
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FIG. 1: DWs for two 2D samples with different bonds distributions 
and periodic BCs in both directions. Points are spins. The lines cross 
the bonds (not shown) that have 'changed state' when antiperiodic 
BCs are introduced in the horizontal direction. The thick line is the 
DW. Left: Gaussian. Right: Bimodal. 
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TABLE I: Parameters used in the MC runs for 2D systems (See text 
for details). 
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TABLE II: Parameters used in the MC runs for 3D systems (See text 
for details). 
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ary conditions (BCs) in all directions, and in 2D we have stud- 
ied two cases: periodic -periodic BCs and periodic-free BCs. 

For the EAG the definition of the DW is straightforward: 
given a system we 'perturb' it by changing the sign of all the 
bonds in a column in 2D or plane in 3D (i.e. we introduce 
antiperiodic BCs in one direction) and compare the GS of the 
new system with that of the unperturbed one. The DW is sim- 
ply the set of bonds which have 'changed state', i.e. those for 
which JijOiOj has changed its sign. It can be shown that in 
the dual lattice, the dual of the bonds that have changed can 
only form a path in 2D (or surface in 3D) 14 of length (or area) 
I, that goes from one border of the lattice to the opposite one 
(see Fig.[T|i. dt is obtained from the relation Ijyw = (I) ~ Lrf, 
where low is the sample average ((• • •)) of /. 

For the EAB one has to be a bit more careful because the GS 
is degenerated. If one performs the same procedure as above, 
but now with a random pair of GSs, one sees that many sets of 
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FIG. 2: (color online) DW length for both Gaussian and bimodal dis- 
tribution of bonds in 2D and 3D for systems with periodic boundary 
conditions in all directions. The exact values were obtained with a 
branch-and-bound algorithm. Error bars are smaller than the sym- 
bols. 



bonds have changed state (see Fig.Q]). Only one of these sets 
forms a structure that crosses the sample, as in the EAG. The 
other sets form loops of zero energy: they enclose clusters of 
spins that can be flipped without changing the energy of the 
sample. As these clusters are present both in the perturbed and 
the unperturbed systems, it is reasonable to define the DW of a 
given pair of GS, as only the structure that crosses the sample. 
In this case the characteristic length of the DW of a sample, 
/, is defined as the average over all pair of GSs of the length 
of DWs, as defined above. In practice, to extract the DW in 
the EAB one has to look for the percolating structure in the 
dual lattice, both for 2D and 3D. For this we have used the 
algorithm of Hoshen-Kopelman. 1 ^ We have been careful to 
eliminate the loops (or closed surfaces) that stick to the DW 
(see Fig.[T]for an example). 

For small EAB samples we have obtained the exact value 
of the DW length by averaging over all pairs of GSs for each 
sample. To obtain all the GSs for systems with fully periodic 
BCs we have used a branch-and-bound algorithm. 16 We have 
studied sizes up to L = 8 in 2D and L = 4 in 3D (10 4 samples 
for each size). Results are shown in Fig. [2] This approach 
is not practical for larger sizes because of the fact that the 
number of GSs grows exponentially with L. 17 Therefore, to 
estimate the average DW length for larger samples, we have 
resorted to a Monte Carlo (MC) algorithm to obtain one pair 
of random GSs for each sample (see Fig.|2|. 

The MC algorithm we have used is a variant of parallel 
tempering 18 but suitably modified to find a ground state as 
quickly as possible. As in Ref. 18, we use a compound sys- 
tem or ensemble, that consists of m noninteracting replicas 
of the system, each one associated to a different temperature 
in the interval [T m j„, T max ], where the distance between con- 
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secutive temperatures is a constant. In general the heart of a 
parallel tempering algorithm consists of two routines that are 
performed alternately. One of them consists of a standard MC 
algorithm applied to each replica: in each elementary step, 
the update of a random selected spin of the ensemble is at- 
tempted with probability given by the Metropolis rule^. In 
the other routine, an exchange of two configurations between 
two replicas at consecutive temperatures is attempted with the 
probability defined in Ref. [HI In general, the unit time or 
MC step (MCS) of a parallel tempering algorithm consists of 
a fixed number of elementary steps of standard MC, followed 
by another fixed number of trials of replica exchange. To equi- 
librate the system a MCS defined as m x N elementary steps 
of standard MC and only one replica exchange, is usually cho- 
sen. As we are only interested in reaching quickly a configura- 
tion of GS, we have used a different MCS: it consists of mxN 
cycles, where a cycle is defined as only one elementary step 
of standard MC plus one replica exchange. After t MCS, the 
algorithm stops and the configuration with minimum energy 
(among all configurations visited in all replicas in the simu- 
lation process) is stored. We have found that this algorithm 
is very efficient for reaching the GS: we have checked that, 
for the sizes studied with the parameters indicated in Tables I 
and II, our heuristic outputs a true GS with a probability larger 
than 0.99. In particular for 2D systems, we have verified this 
by calculating the ground states energies with an exact branch- 
and-cut algorith m 2Q i 21 and comparing them with the energies 
of the configurations we obtained (a detailed analysis of this 
algorithm will be published elsewhere—). 

In our simulations, for each lattice size L of EAB and EAG 
in both 2D and 3D systems we have used m — 20 replicas with 
temperatures between T max = 1 .6 to T m i„ = 0.1. The number 
of MCS t and the number of samples Ns analyzed for each 
sample size, is given in Tables I and EL 

Interestingly, considering only one pair of GSs for each 
EAB sample leads to an estimate of the average DW length 
that is within only 1% of its exact value (for small systems). 
This is good evidence that the MC algorithm is sampling the 
GS space almost uniformly. As our algorithm does not reach 
equilibrium we have checked explicitly that the sampling is 
indeed almost uniform. In Fig. [3] we show the distribution of 
GSs reached by our algorithm for a typical 2D sample with 
80 GSs. It is instructive to compare our method with the sam- 
pling that results from the algorithm (hereafter called CBA) 
used by Cieplak and Banavar. 23 In this technique an infinites- 
imal noise is added to the couplings, and the GS of the system 
is obtained. Then, the same realization of the noise is added 
to the system with antiperiodic BCs and a new GS is obtained. 
The noise in the couplings breaks the degeneracy, and the fact 
that the same realization of the noise is used ensures that the 
resulting DW has no loops. Up to L = 8 we have also used 
the CBA to analyze the same samples as with the branch-and- 
bound algorithm. As Fig.|3]shows, the sampling given by this 
technique is far from uniform. In fact, there is a large fraction 
of GSs that are not reached by it. In spite of this, the estimate 
of the average DW length given by the CBA is within 3% of 
its exact value (it must be stressed that the bias is systematic: 
the estimates are consistently smaller than the exact values). 
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FIG. 3: (color online) Histogram of the fraction of runs for which 
each GS has been obtained, for two different stochastic algorithms 
(see text). The states have been ordered by decreasing frequency 
of sampling, for each algorithm, i.e. the orderings depend on the 
algorithm. In the inset the ordering is the same for both algorithms. 



This suggests that the sampling method might not be crucial 
to obtain good estimates of df. 

For 3D EAG systems (see Fig. |2]i the value we obtain is 
df = 2.68(2). This value is compatible with the ones found by 
other authors: df = 2.68(2) and df « 2.7 i 24 ' 25 However, the 
fractal dimension obtained when comparing GSs obtaining by 
studying the response of the GS to a coupling-dependent bulk 
perturbation is a bit smaller: df — 2. 57(2). 26 This can perhaps 
be explained by arguing that, given the fact that droplets of 
all sizes appear when there is a bulk perturbation, the fractal 
dimension of all the droplets do not have to coincide. Recently 
it has been found that in 2D droplets with different sizes do 
seem to have different values of df. 21 

For 3D EAB systems we have compared the average DW 
lengths obtained using our MC algorithm and the CBA algo- 
rithm, up to L = 7. As in the 2D case, the difference between 
both values is almost constant and smaller than 1 %. Using our 
MC algorithm up to L = 10 we obtain df = 2.70(2). A differ- 
ent value is obtained in Ref. [TH where the DW is restricted to 
the rigid lattice. This is the set of bonds which do not change 
its state in all the GSs (i.e. the bonds are either frustrated 
or unfrustrated in all GSs). The value reported there for 3D 
(df = 2.59(2)) is a bit smaller than the one reported in this 
article. This is due to finite size effects, given the fact that the 
size of the rigid lattice varies greatly from sample to sample. 

For 2D systems we obtain df = 1.33(1) for the EAB and 
df = 1.27(1) for the EAG (see Fig.|2|. Notice that these two 
values are clearly different. Given the small systems consid- 
ered, this discrepancy could be attributed to finite size effects 
(which are usually large in 2D systems). As much larger sizes 
can be analyzed in lattices with one free BC, we have stud- 
ied the length of the DW in those systems. It is well known 28 
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that in this case the problem of finding the GS can be mapped 
to a minimum-weighted perfect matching problem, for which 
very efficient algorithms exist. We have used one implementa- 
tion of a Blossom algorithm, 29 which has allowed us to obtain 
the GSs up to L = 300 for EAB and up to L = 100 for EAG 
(10 3 — 10 4 samples for each size). 

As in the case of systems with fully periodic BCs, for the 
EAB we have calculated the exact average length of the DW 
for small systems (up to L — 16). For this we have used an al- 
gorithm to count all the minimum- weighted perfect matchings 
of each sample. 30 We have also applied the MC algorithm, 
which has allowed us to reach sizes up to L = 22. Again, the 
MC estimate is within 1% of the exact value. To reach even 
larger sizes we have used an algorithm based in the Blossom 
routine to choose a pair of random GSs. This was done by 
choosing a GS among all the GSs compatible with a given ran- 
dom optimal matching. The results are shown in Fig. [4] In the 
inset the three approaches are compared, for small systems. It 
can be seen that the estimate given by the Blossom algorithm 
is clearly smaller than the exact value. This is evidence that 
not all GSs are reached with the same probability. 

To understand the origin of this bias it is necessary to no- 
tice that each optimal matching corresponds to many GSs. We 
call this the degeneracy of the matching. But these degenera- 
cies vary widely from matching to matching. Therefore, a 
random sampling of the optimal matchings tends to favor GSs 
present in matchings of small degeneracies. For some reason, 
these particular GSs of the periodic and antiperiodic system 
are more similar than GSs chosen uniformly at random. This 
leads naturally to smaller DWs. However, the bias introduced 
by this nonuniform sampling is so small (less than 6%) that it 
is difficult to make this reasoning more precise (see the inset 
in Fig. |4). In spite of this bias, the exponent estimated with 
this procedure is very close to the value given by the exact 
and MC algorithms, up to L = 22. Thus, as in the case of our 
MC algorithm and the CBA, the only effect of the bias is to 
shift the points, but without changing the slope of the fit. This 
shows again that the sampling method does not seem to be 
crucial to determine the fractal dimension. 

Studying systems of sizes up to L = 300 the value that we 
estimate for the fractal dimension of the 2D EAB model is 
df = 1.323(3). This value is between the bounds given in 
Ref. and is compatible both with the value obtained for the 
case of fully periodic BCs, and with the most recent estimate 
calculated using the entropy ansatz: df = 1.30(3)^ Even 
though a value of df — 1.30(1) has been reported previously, 12 
the small discrepancy with the value reported in this letter is 
due to the larger system sizes that are considered here (up to 
L = 100 in Ref.[l2r). For the sake of comparison, we have also 
calculated the fractal dimension for the 2D EAG. We obtain 
df = 1.279(4), which is compatible with the most recent and 
accurate value reported: df = 1.274(2). 7 

Notice that the fractal dimensions obtained for the 2D EAB 
and the 2D EAG are clearly different. To check that this is not 
an artifact of the scaling function used to fit the points, BL /, 
we have also tried fits with two other scaling functions with 
additional correction terms (see Table III). For EAB we have 
obtained df > 1.32 for all the fits with Q > 0.4 and Q > 0.1 
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FIG. 4: (color online) DW lengths for 2D systems with periodic-free 
BCs. The dashed line represent the bounds found in Ref. 7. The inset 
shows both the exact average and the Blossom estimate of the DW 
length, as well as a Monte Carlo estimate, for samples of small size. 
Error bars are smaller than the symbols. 



for periodic-free and periodic -periodic systems, respectively 
(a fit with Q > 0.1 is considered as a good fit 31 ). This shows 
that, with a high probability, the difference found between the 
fractal dimensions for the EAB and the EAG is not an effect 
of finite size scaling. This difference is further evidence that 
the universality of both models is different for T = 0. 

As mentioned above, Eq. ((TJ (considering 0=0) predicts 
jSLE _ 1 25. This is the same universality class as the loop- 
erased random walks. 32 Interestingly, the fractal dimension 
we obtain is much closer to df =4/3, which is the fractal 
dimension of the normal self-avoiding walks. 33 This does not 
necessarily mean that DWs cannot be described as SLE pro- 
cesses. The walks generated by these processes have many 
interesting properties. One way to see whether DWs can be 
described as SLE processes is to test whether they have some 
of these properties. 5 

We have performed one such test in 2D. If DWs are de- 
scribed by an SLE process of diffusion constant K, the proba- 
bility that they pass to the right of a point with polar angle 
(measured from the starting point of the walk) is- 3 - 4 -: 



i r(i)cot(0) 



1 4 3 
2' Jc'2 



;-cot 2 (0) 



(2) 



where 2-F1 is the hypergeometric function. Notice that Pk(<I>) 
does not depend on the radial coordinate R of the reference 
point. The result of this test is shown in Fig. [5] The points for 
the EAG are compatible with a diffusion constant k = 2.23, 
which is consistent with the value obtained by Bernard et al. A 
In turn, this value is compatible with the relation df = 1 + if/8. 
For the EAB the points obtained are far from the curve (upper 
curve in Fig. [5]) that corresponds to K — 8(df — 1) = 2.584 
(using df = 1.323). In fact, they are much closer to ^2.23(0) 
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TABLE III: Fractal dimension obtained using different scaling corrections. The fits are restricted to the interval [L,„i n ,L max ]. Q gives the 
quality of the fits (see the text). 
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FIG. 5: (color online) Comparison between the prediction of Eq. l(2j 
and the values obtained numerically for 2D EAB and 2D EAG sys- 
tems. The lines show the difference between Pk{<P) and /*>(0) for 
two values of K. The points are averages over 5000 samples for EAB 
systems with L = 100 and R = 30, and over 2000 samples in the other 
cases. 

(lower curve in Fig. [5). Thus, unless finite size effects are 
very significant, this result shows that EAB DWs cannot be 
consistently described as SLE processes. It must be stressed 
that this test has been performed using the DWs calculated 
with the Blossom algorithm which means that, for the EAB, 
the GSs have not been sampled uniformly. Nevertheless, as 



the determination of the fractal dimension of the DWs does 
not seem to be affected by this bias, it is not unreasonable to 
assume that it also does not affect the outcome of the SLE 
test. We have only used the Blossom algorithm because, for 
this kind of tests, very large system sizes must be studied to 
obtain meaningful results. 

In conclusion, we have performed an extensive numerical 
study of the domain walls for the EAB and the EAG both 
in 2D and 3D, and have determined their fractal dimension 
from a direct measurement of their length. In 3D we find 
that both exponents coincide. On the other hand, for 2D sys- 
tems we find a clear difference between the fractal dimen- 
sions of the EAB and the EAG. The value obtained for the 
EAB model shows that these domain walls are more similar to 
self-avoiding walks than to loop-erased random walks, as one 
would expect from Eq. (|T). Even though the methods we have 
used do not sample uniformly the ground state space, the fact 
that all the obtained values of df coincide is good evidence 
that this estimate must be very close to its true value: given 
the very different nature of the algorithms used it would be 
very unlikely that they give raise to the same wrong estimate 
of the fractal dimension. 35 Finally we have tested whether the 
2D EAB domain walls can be described by SLE processes. 
The outcome of this test is negative: the probability that a 
DW passes to the right of a given point is not consistent with 
an SLE process with K = %{df— 1). 
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